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0> Abstract 

m 

In this paper, we develop modified versions of the likelihood ratio test for multivariate het- 
eroskedastic errors-in-variables regression models. The error terms are allowed to follow a mul- 
tivariate distribution in the elliptical class of distributions, which has the normal distribution as 
a special case. We derive the Skovgaard adjusted likelihood ratio statistics, which follow a chi- 
squared distribution with a high degree of accuracy. We conduct a simulation study and show that 
the proposed tests display superior finite sample behavior as compared to the standard likelihood 
ratio test. We illustrate the usefulness of our results in applied settings using a data set from the 
WHO MONICA Project on cardiovascular disease. 
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1 Introduction 



In actual problems, it is not uncommon to observe variables subject to measurement errors. This 
characteristic leads to extra variations not predicted by the standard model. Such extra variations 
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must be taken into account in statistical modeling to avoid misleading conclusions (e.g., inconsistent 
estimators, false confidence intervals, and underestimated standard errors). Fuller [13J presents 
many interesting problems with this peculiarity, one of his examples is presented in the following. 
Suppose that we wish to fit a regression line for corn yield against nitrogen in the soil. It is known 
that the true quantity of nitrogen in the soil is not directly available and it has to be estimated 
by laboratory analysis from a selected sample. In order to capture this extra variability of the 
collected data, measurement error equations are included in the modeling process. This procedure 
produces, in general, consistent estimators for the model parameters. Therefore, the full specification 
of errors-in-variables regression models is essentially based in two sets of equations, one related to 
the regression modeling (known as the regression equation) and other for the measurement errors 
(known as measurement equation). 

Measurement error models are widely applied in many fields, such as epidemiology, astronomy, 
biology, and so on. In the following we present some applications that motivated us the development 
of the present article. Aoki et al. p] employ a statistical analysis to compare the effectiveness of 
two types of toothbrushes in removing plaque; the covariate (independent variable) is a plaque score 
before brushing and the response variable (dependent variable) is the plaque score after brushing. 
The amount of plaque before and after brushing is obtained in the same way. In this case, as the 
amount of plaque is imprecisely evaluated, it is reasonable to consider that the covariate is subject to 
measurement errors. Kulathinal et al. |15j and Patriota et al. |20j study a univariate heteroskedastic 
errors-in-variables model with equation error and normal distribution. An extension to the scale- 
mixture of normal distributions is presented by Cao et al. [5j. The results in these three latter 
papers were applied to an epidemiological dataset of the WHO MONICA (World Health Organization 
Multinational MONItoring of trends and determinants in CArdiovascular disease) Project. Patriota 
et al. |21j extend this univariate version to the multivariate one keeping the normality assumption, 
and also develop bias correction for the maximum likelihood estimators. 

In the aforementioned problems, the appropriate model should consider the measurement equation 
in the modeling, as presented in Section [2j When the true covariates (the unobservable covariates) 
are modeled as random variables, the model is said to be structural. Otherwise, when these true 
covariates are fixed unknown parameters, the model is said to be functional. The structural model 
has identifiability problems while the functional version has unbounded likelihood. To overcome such 
problems, it is usual to regard some quantities as known, e.g., the measurement errors variances ([7], 
[25]) or the model intercept pQ. Here, we adopt the structural version; for more details regarding 
errors-in-variables models we refer the reader to Fuller [13] , Cheng and Van Ness [8j and Buonaccorsi 

Generally speaking, inference procedures under complicated models are based on asymptotic 
arguments. However, as shown by Melo and Ferrari [18], first-order asymptotics under errors-in- 
variables models may be inaccurate for small- and moderate-sized samples. With the purpose of 
improving the approximation of the distribution of the likelihood ratio statistic to its asymptotic 
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reference chi-square distribution, Skovgaard [23] proposes a new test statistic when the interest lies 
in testing a parameter vector (see Section 3 for details). In this paper, we apply his proposal to 
a multivariate regression model with measurement errors when the random quantities have a joint 
elliptical distribution. The model considered here is a generalization of the one studied by Patriota 
et al. [21]. Therefore, our results apply to this latter paper and to all its subcases. 

It is well-known that the normality assumption is not always tenable and alternative distribu- 
tions (or methodologies) should be considered in such situations. For this reason, we consider the 
elliptical family of distributions, since this class includes many important distributions such as nor- 
mal, Student-t, power exponential, contaminated normal, Pearson II, Pearson VII, and logistic, with 
heavier or lighter tails than the normal distribution. Basically, we say that a q x 1 random vector 
Z has a multivariate elliptical distribution with location parameter /x (q x 1) and a positive definite 
scale matrix fi (q x q) if its density function exists, and is given by Fang et al. 

f z (z) = \n\- 1 i 2 Po [(z- n) T or\z-p)\, 2eR«, (l) 

where po : E — > [0, oo) is such that J °° uz~ 1 po(u)du < oo. We use the notation Z ~ El q (fi, fi,po). 
It is possible to show that the characteristic function of Z is given by if>(t) = E(exp(it T Z)) = 
exp(ii T /x)0(t T fit), where t £ R' and <ft : [0, oo) — > M. If 4> is twice differentiable at zero, then 
E{Z) = n and Va,i(Z) = cfi, with c = —2cp'(0). For more details, we refer the reader to Fang et al. 
[TT] and Fang and Anderson [12]. 

The paper unfolds as follows. Section [2] introduces the heteroskedastic multivariate regression 
model with measurement error. Section [3] contains our main results, namely the ancillary statistic 
and an explicit formula for the modified likelihood ratio test. Section [4] presents a simulation study 
on the finite sample behavior of the standard likelihood ratio test and its modified counterpart. Our 
simulation results show that the likelihood ratio test tends to be liberal and its modified version is 
much less size-distorted. An application that uses real data is presented and discussed in Section [5] 
Finally, Section [6] concludes the paper. Technical details are collected in three appendices. 

2 Heteroskedastic multivariate regression model with measurement 
error 

We consider the following heteroskedastic multivariate regression model with measurement error 

yi = a + (3xi + qi, 

Yi = yi + a, (2) 

for i = 1,2, ... ,n, where yi is an (m x 1) latent response vector, Xi is a (p x 1) latent vector of 
covariates, a is a vector of intercepts, and /3 is an (m x p) matrix of parameters. The error terms 
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are qi,ei and itj. The vectors yi, and X( are not directly observed; instead, the vectors Yi and 
X{, respectively, are observed. The observed and unobserved vectors are connected by the additive 
relation specified by the two last measurement equations in (12]). Model Q can be shortly written as 



Zi = 



S + Abi 



(3) 



with % = 1, 2, . . . , n, where 



Zi = 



( 




/3 Im Im 

Ip o I p 



) 



and bi 



Bi 



\UiJ 



where I p is the identity matrix of dimension p x p. We assume that b\, . . . , b n are independent 
and identically distributed according to an elliptical distribution with the function <p being twice 
differentiate at zero (see the previous section, below Eq. Q). This model specification generalizes 
the one proposed by Patriota et al. |21j . Here, we define the covariance matrix of bi as 



do not depend on the index i then we have the homoskedastic multivariate regression model with 
measurement error. In short bi ~ El^ 2m +2p){ T 1^bi'iPo)i with 



and c = —2(p'(0). It can be shown that the vector Zi has a (p + m)-variate elliptical distribution with 
location vector fi = £ W +m and dispersion matrix $7j = Qi(6). The parameter vector is 6 = 
(vec(/3) T ,Q; T ,/iJ,vech(Sg) T , vech(Sa;) T ) T £ 6 C I s , where s = mp+m+p+p(p+l) /2 + m(m + l) /2 
is the dimension of the parameter space, vec is the vector operator stacking the columns of a matrix 
underneath each other, and vech is the operator that stacks the lower triangle of a matrix into a 
column vector. By using properties of elliptical distributions, we can write Zi ~ El p+m (fi,£li;po), 
where /i = 5 + Ar] and Qi = AO{,. A T are the location vector and scale matrix, respectively, which 
are functions of the parameter vector 6, namely: 



/£* \ 




77 = 








n bt = c - l iii 




and 
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Notice that, for the normal case, model ([3]) is identifiable. Hence, if the density generator po does 
not depend on additional parameters to be estimated, model ([3]) is also identifiable for the class of 
elliptical distributions. It is noteworthy that for many elliptical distributions po contains unknown 
quantities, e.g., the degrees of freedom in the Students distribution and the shape parameter in the 
power exponential distribution. We may want to estimate these quantities via maximum-likelihood 
estimation, however, Lucas [T7] studies some robust features of M-estimators by using influence 
functions for the Student-i distribution and shows that the protection against "large" observations is 
only valid when the degrees of freedom parameter is known. Therefore, for the purpose of avoiding 
possible lack of protection against outliers, we do not estimate those unknown parameters. Otherwise, 
unboundedness problems may arise for the influence functions and the elliptical distribution will lose 
its main goal. In practice, we can use model selection procedures, such as the Akaike information 
criterion (AIC), to choose the more appropriate values of such unknown parameters. 

The likelihood function for model O) is 



Lemonte and Patriota [16] discuss the maximum likelihood estimator of 9, and derive the gener- 
alized leverage as well as the normal curvatures of local influence under some perturbation schemes. 
The reader is referred to their work for further details on the iterative procedure to attain the 
maximum likelihood estimates. 

3 Modified likelihood ratio test 

Let 6 = (ip T , lj t ) T , where ij) is the vector of parameters of interest and u; is the vector of nuisance 
parameters. Here, the null and alternative hypotheses are, respectively: Hq : ip = tp^ and Hi : ift ^ 
where is a known q- vector. In this paper, we focus on testing hypotheses on some elements 
of the matrix f3. Define vec(/3) = (fil , ^J)" 1 "; without loss of generality, we consider tp = (3± and 
u = (f3j , ct T , /xj, vech(S g ) T , vech(S x ) T ) T . The maximum likelihood estimator over O is denoted by 
6 = (tp,LJ T ) T and the maximum likelihood estimator under the null hypothesis, by = (t/>,c5 T ) T , 
where ip = We use and "~" for matrices and vectors to indicate that they are computed 

at 6 and 9, respectively. 



n 



L(e)= Pz (z,e) = l[\a 



where di = di(6) = Zi — \x. The logarithm of the likelihood function is 




(4) 
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The likelihood ratio statistic for testing T~Lq is 



LR 



£$)-£$)}. 



Under regularity conditions on po see, for instance, Cox and Hinkley [9], this statistic converges in 
distribution, under Hq, to X% , where q is the dimension of ip. 

When the parameter of interest, if), is one- dimensional, we can use a signed likelihood ratio 
statistic to test Hq against %\. Barndorff- Nielsen [3] proposes a modified version of this statistic 
that depends on an ancillary statistic and some derivatives with respect to the sample space. This 
modified version intends to better approximate the signed likelihood ratio statistic distribution for the 
standard normal distribution. The extension for the vectorial case was developed by Skovgaard [23], 
for which two modified statistics, that are asymptotically equivalent, were proposed. We denote these 
two statistics here by LR* a and LR*,* . These modified statistics retain the essential character of LR, 
but can be difficult to obtain. The difficulty arises from the need of an appropriate ancillary statistic 
and derivatives of the log-likelihood function with respect to the data. By "ancillary statistic" we 
mean a statistic, say a, whose distribution does not depend on the unknown parameter 0, and that 
along with the maximum likelihood estimator 0, is a minimal sufficient statistic for the model. If 
(6, a) is sufficient, but not minimal sufficient, Barndorff-Nielsen's results still hold; see, Severini [22] 
(§ 6.5). In fact, minimal sufficiency is only required for the ancillary a to be relevant to the statistical 
analysis. Sufficiency implies that the log-likelihood function depends on the data only through (6, a), 
and we then write 1(0; 6, a). The required derivatives of £(0; 9, a) with respect to the data are 

/ = 31(0; 0, a) u> = d 2 e(0;d,a) ^ ~ _ d 2 £(G;G,a) 
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The modified likelihood ratio statistics are 

1 



0=0, 0=0 



LR* a = LR[l-—\ogp 
and 

LRl* = Li? -2 log p, 

with 

p = | J |1/2|^/|-1| j |1/2|J |-l/2|j |l/2 ! ^ (5) 

LRi/ 2 - 1 (e' - l') T (U')- l U 

The quantities £' — £', U' and J = J(0;0,a) are derivatives with respect to the sample space and 
are obtained from an ancillary statistic a and the maximum likelihood estimator 9. Matrix 
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is the lower-right sub-matrix of J related to the nuisance parameters u>. Although the statistic 
LR* a is non- negative, and when q = 1, it reduces to the square of Barndorff-Nielsen's statistic, the 
second version, LR** , seems to be numerically more stable and is naturally attained from theoretical 
developments see Skovgaard [23]. These statistics can be approximated to with high accuracy 
under the null hypothesis. 

In order to find an ancillary statistic, we first note that model ^ is a transformation model. Hence 
any maximal invariant statistic is an ancillary statistic; see Barndorff-Nielsen et al [1] (Chapter 8) 
and Barndorff-Nielsen [2j. Moreover, a statistic is maximal invariant if it is a function of any invariant 
statistic. Pace and Salvan jTU] (Theorem 7.2) showed that in transformation models all invariant 
statistics are distribution constant, i.e., their distributions do not depend on 6. Let Pj = Pi(9) be a 
lower triangular matrix such that P%Pj = is the Cholesky decomposition for alH = 1, . . . , n. By 
using the same idea as Melo and Ferrari [18], the statistic a = (aj , aj, . . . , a,[) T , with 

a; = Pr 1 ( Zi - fi) , (6) 

is maximal invariant and consequently is ancillary, as previously shown, where p, = /x(0) and 
Pi = Pi{6) are the maximum likelihood estimates of /x and Pi, respectively. Taking this ancillary 
statistic, we can find the derivatives with respect to the sample space required for Skovgaard's 
adjustment. 

In order to compute the score vector, observed Fisher information, and the derivatives with 
respect to the sample space, we define the following matrices and vectors. Let T be a matrix, with 
its (j,k)-th element given as 

n 

where fig. = d£li/d8j and fie-e fc = d£lgJd6}-, for j,k = 1,2, ...,s. Observe that fig. does not 
depend on the index i, since the matrix fij merely depends on i through the known matrices S ei , 
£( ue v, and S Uj . Therefore, their derivatives with respect to the model parameters are zero. Consider 
the following block- diagonal matrices R = diag(r, r, . . . , r) and V = diag(t>, v, . . . , v), of dimension 
sn x s. The i-th elements of r and v are, respectively, r$ = W Po (dJ di) and v% = W'(djQ,~ di). 

Let h = (^h^ T , . . . , h^ T ^j and w = (w^ T , . . . , w^ T ^j be sn-vectors, for which the i-th elements 

are given, respectively, by = — dTfi-ifi^.fir 1 ^ - 2/xJ.fir 1 d l , and wf^ = (Pe^i + fie^j ^ 1 d i , 
where Pg. = dPi/d9j and fig. = dpi/dOj. In addition, let B, C, F, G, M, and Q be (sn x s) 
block-matrices, whose (j, k)-th blocks are, respectively, given by the vectors b^ k \ c^ k \ f^ k \ g^ k \ 
m^ k \ and q^ k \ The i-th components of these vectors are, respectively, given by b^ k ' = h^'h^, 

c i = w i h i > 
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mp' fc) = dj (2nr 1 n 6k n- i n g .nr 1 - n^n^nr 1 ) d % + 2^ or^.or 1 ^ + 2/xJ.nr 1 n 9k nr 1 d i 

where /^e^ = d^iQ .jdQ^. The score vector, observed information matrix, and derivatives with respect 
to the sample space are given in matrix notation by 

U = -- n* + R T h, J =^T - R T M -V T Q, 

i' = 2R T w, U' = 2{R Y B + V T C^j , and J = 2 (r t F + V J G) , 

where the j-th component of the vector n* is Y^l=i * r (^~ ^%) - 

When we insert J, J W a>j J^w, U', £' — £', and the original likelihood ratio statistic LR into 

([5]), we obtain the required quantity p for Skovgaard's adjustment. Now, one is able to compute the 
modified versions LR* and LR** . 

Computer packages that perform simple operations on matrices and vectors can be used to cal- 
culate p. Note that p depends on the model through /x, Pj, Oj and fij" . The dependence on the 
specific distribution of z in the class of elliptical distributions occurs through W Po . Appendix A gives 
the required derivatives of \x and f2j with respect to the vector of parameters. The derivative of P% 
with respect to the vector of parameters was obtained by using the algorithm proposed by Smith 

4 Simulation study 

In this section, we present simulation results to illustrate the effectiveness of the Skovgaard's adjust- 
ments for the original likelihood ratio statistic. The performance of the statistics LR, LR* a , and LR** 
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were evaluated considering the repeated sampling principle. We study the frequencies of simulated 
samples that commit the type I error, and compare them with actual pre- assigned nominal levels. 

The simulation study was based on model ^ when m = 1 and the random (p + l)-vector Z{ 
follows a multivariate normal, Student-t with v degrees of freedom or power exponential distribution 
with shape parameter A (Gomez et al. [14J). This study was conducted using the matrix programming 
language Ox Doornik [10J considering 10,000 (ten thousand) Monte Carlo samples. We consider the 
following sample sizes: n = 20,30,40 and 50. The nominal levels are 7 = 1%, 5%, and 10%. In 
addition, we consider two to four covariates, i.e., p = 2, 3, and 4. The null and alternative hypotheses 
are, respectively, Hq : tj) = and %\ : -0 ^ 0, where rj) = (/3i,/02> • • • > Aj) T , for q = 2,3 or 4. The 
true parameters are a = 0.2, fi x = (—2, —2, . . . , — 2) T , T, q = 10 and T, x = 4I p , where I p is the (p x p) 
identity matrix. As before S ej , £( ue v and for i = 1, 2, . . . , n, are assumed to be known. We set 
S( ue v as a p x 1 vector of zeros, the values for T, ei , for i = 1, . . . , n were generated as n independent 
observations from y^S^~ ~ U(0, 1), and T, Ui , for i = 1, . . . n, was assumed to be a diagonal matrix; 
the square root of its elements were independently generated from a U(0, 1) distribution. After 
generating these scale matrices, they were kept fixed for all Monte Carlo Simulations. For the 
Student-i distribution, we fixed the degrees of freedom at v = 5 and, for power exponential model, 
the shape parameter was fixed at A = 0.6. Tables 1, 2 and 3 present the rejection rates for the tests 
based on LR, LR*, and LR^* , for q = 2 and different values for p and n. 

In our simulations we observe that tests based on the original likelihood ratio statistic LR are 
liberal when the sample size is small, since the respective rejection rates are larger than the actual 
nominal levels. This can be seen in Table 2 when Zi follows a Student-i distribution, p = 3, q = 2, 
and n = 20. In this case, the rejection rates for tests based on LR equal 2.5% (7 = 1%), 9.0% 
(7 = 5%), and 15.7% (7 = 10%). When Z; L follows a normal distribution, p = 2 = q, and n = 20, 
we observe the following rejection rates: 2.1% (7 = 1%), 8.2% (7 = 5%), and 14.7% (7 = 10%); see 
Table 1. 

On the other hand, tests based on the modified versions LR* a and LR^* present rejection rates 
closer to the actual nominal levels than the original version LR. Table 2 gives the rejection rates 
when Zi follows a Student-t distribution, p = 4, and q = 2. For n = 20 and 7 = 10%, the rejection 
rates are 16.1% (LR), 10.9% (LR* a ), and 10.2% {LR**). For n = 30 and 7 = 10%, the rejection rates 
are 13.9% (LR), 10.2% (LR* a ), and 10.0% (LR**). Another example can be seen in Table 1 for the 
normal distribution, p = 3, q = 2, n = 30, and 7 = 1%. For this case, the rejection rates are 2.2% 
(LR), 1.2% (LR*,), and 1.2% (Li?**). For the power exponential distribution, p = 4, q = 2, n = 30, 
and 7 = 5% we have 7.8% (LR), 5.1% (LR*), and 4.9% (LR**); see Table 3. 

Figures [TJ [|and [3] depict curves of quantile relative discrepancies versus the correspondent asymp- 
totic quantiles for the statistics LR, LR*,, and LR*J* considering different sample sizes and p = 4 
and q = 3 under a Student-t distribution (Figure [l| , p = q = 4 under a normal distribution (Figure 
[2]), and p = q = 3 under a power exponential distribution (Figure [3]). The closer to zero a curve is, 
the better is the asymptotic approximation. As expected, the three figures show that the curves of 
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relative quantile discrepancies for the statistics LR* a and LR*,* are closer to the horizontal axis of the 
zero ordinate than the curve based on LR. 

Tables 4 and 5 present rejection rates obtained under alternative hypotheses for n = 20 and 
n = 40, respectively. Recall that the null hypothesis is "Ho ■ ip = 0, and that tp = (/3±, . . . , (3 q ). Here, 
the rejection rates were obtained assuming that i/> = (77, . . . , 77), with r\ varying from 0.1 to 1.5, and 
p = q = 2. The original likelihood ratio test was not included in this study, since our simulations 
indicated that it is oversized. The figures in Tables 4 and 5 indicate that the modified tests present 
similar power, the power of the test that uses LR* a being equal to or slightly greater than that of the 
power of the test based on LR** . 

We conclude that the modified versions of the likelihood ratio test perform better than the original 
test for small and moderate sample sizes. Although the tests based on the modified statistics, LR*, 
and LR*,*, present similar results, the test based on LR*,* has a slightly better performance (in the 
majority of the cases) than the one based on LR*,. 

Table 1: Rejection rates of when q = 2 and p = 2,3, 4; normal distribution. 



p = 2 






7 = 1% 






7 = 5' 


Vo 




7 = 10% 


n 


LR 


LRl 


LR** 


LR 


LR* a 


LR* a * 


LR 


LRl LR* a * 


20 


2.1 


1.1 


1.0 


8.2 


5.4 


5.2 


14.7 


10.6 10.4 


30 


1.7 


1.1 


1.1 


6.9 


5.2 


5.1 


13.0 


10.6 10.5 


40 


1.5 


1.1 


1.1 


6.2 


4.9 


4.9 


11.9 


10.0 10.0 


50 


1.5 


1.1 


1.1 


6.2 


5.2 


5.2 


11.7 


10.3 10.3 


p = 3 






7 = 1% 






7 = 5' 


Vo 




7 = 10% 


n 


LR 


LRl 


LR* a * 


LR 


LRl 


LR* a * 


LR 


LRl LR* a * 


20 


2.6 


1.1 


1.0 


9.3 


5.6 


5.3 


16.0 


10.5 10.2 


30 


2.2 


1.2 


1.2 


7.8 


5.3 


5.1 


13.9 


10.6 10.4 


40 


1.5 


1.0 


1.0 


6.9 


5.1 


5.0 


12.8 


10.4 10.3 


50 


1.6 


1.1 


1.1 


6.8 


5.4 


5.4 


12.1 


10.4 10.4 


p = 4 






7 = 1% 






7 = 5' 


Vo 




7 = 10% 


n 


LR 


LRl 


LR** 


LR 


LRl 


LR* a * 


LR 


LRl LR* a * 


20 


3.4 


1.2 


1.1 


11.3 


6.1 


5.7 


18.8 


11.7 11.2 


30 


2.1 


1.1 


1.0 


8.4 


5.4 


5.2 


15.0 


10.6 10.4 


40 


1.6 


1.1 


1.0 


7.3 


5.1 


5.1 


13.9 


10.7 10.5 


50 


1.5 


0.9 


0.9 


6.5 


4.9 


4.8 


12.7 


10.3 10.3 
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2: Rejection rates of Ho when q = 2 and p = 2,3, 4; Student-i distribution (y 



p = 2 






7 = 1% 






7 = 5' 


% 




7 = 10% 


n 


LR 


LRl 


LRl* 


LR 


LR* a 


LiC 


LR 




20 


1.7 


0.9 


0.8 


7.4 


4.8 


4.6 


13.6 


10.1 9.8 


30 


1.4 


0.9 


0.9 


6.5 


4.7 


4.6 


12.1 


9.9 9.8 


40 


1.4 


1.0 


1.0 


6.1 


4.9 


4.9 


11.5 


9.7 9.6 


50 


1.3 


1.0 


1.0 


6.2 


5.1 


5.1 


11.5 


10.2 10.2 


p = 3 






7 = 1% 






7 = 5 


% 




7 = 10% 


n 


Li? 


LRl 


LR* a * 


LR 


LR* a 


LR* a * 


LR 


LRl LR* a * 


20 


2.5 


1.2 


1.1 


9.0 


5.5 


5.3 


15.7 


10.9 10.5 


30 


2.1 


1.2 


1.1 


7.4 


5.4 


5.3 


13.5 


10.4 10.3 


40 


1.5 


1.0 


0.9 


6.3 


4.8 


4.8 


12.0 


9.8 9.7 


50 


1.4 


1.0 


1.0 


6.1 


4.8 


4.8 


11.8 


10.1 10.0 


p = 4 






7 = 1% 






7 = 5 


% 




7 = 10% 


n 


Li? 


l^: 


Li?r 


LR 


LR* a 


LR* a * 


LR 


LR* a LR* a * 


20 


2.6 


1.3 


1.0 


9.2 


5.6 


5.1 


16.1 


10.9 10.2 


30 


2.0 


0.9 


0.9 


7.5 


5.2 


5.1 


13.9 


10.2 10.0 


40 


1.9 


1.2 


1.1 


7.2 


5.5 


5.5 


13.1 


10.5 10.4 


50 


1.5 


1.0 


1.0 


6.6 


5.1 


5.1 


12.3 


10.2 10.1 
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Table 3: Rejection rates of Ho when q = 2 and p = 2,3, 4; power exponential distribution (A 



p = 2 






7 = 1% 






7 = 5' 


Vo 




7 = 10% 


n 


LR 


LR* a 


Li?r 


LR 


LR*. 


LR* a * 


LR 


T E?* T E?** 

LK a LK a 


20 


1.9 


1.0 


1.0 


7.2 


5.0 


4.9 


13.5 


9.8 9.6 


30 


1.4 


1.0 


0.9 


6.4 


4.7 


4.7 


12.2 


9.9 9.9 


40 


1.3 


1.0 


1.0 


5.9 


4.9 


4.9 


11.5 


10.0 9.9 


50 


1.4 


1.2 


1.2 


6.2 


5.1 


5.1 


11.8 


10.4 10.4 


p = 3 






7 = 1% 






7 = 5' 


Vo 




7 = 10% 


n 


Li? 


LR *a 


Li?r 


LR 


LR* a 


LR** 


LR 


LRl LR** 


20 


2.4 


1.2 


i.i 


9.3 


5.5 


5.3 


15.5 


10.9 10.6 


30 


1.9 


1.1 


i.i 


7.5 


5.4 


5.3 


13.1 


10.1 9.9 


40 


1.8 


1.1 


i.i 


6.7 


4.9 


4.9 


12.4 


10.3 10.2 


50 


1.5 


1.2 


1.2 


6.2 


4.9 


4.9 


12.5 


10.3 10.3 


p = 4 






7 = 1% 






7 = 5' 


Vo 




7 = 10% 


n 


Li? 


LR *a 


LR** 


LR 


LRl 


LR** 


LR 


LRl LR** 


20 


2.9 


1.2 


1.0 


10.3 


5.9 


5.4 


17.3 


11.6 10.9 


30 


2.1 


1.2 


1.2 


7.8 


5.1 


4.9 


14.0 


10.2 10.0 


40 


1.9 


1.1 


1.1 


6.9 


4.8 


4.8 


12.8 


10.0 9.8 


50 


1.5 


1.1 


1.0 


6.6 


5.1 


5.1 


11.9 


9.8 9.7 
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n = 20 
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igure 1: Quantile relative discrepancies: p = 4, q = 3, and Student-i distribution for Z{ 
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Figure 2: Quantile relative discrepancies: p = 4, q = 4, and normal distribution for Z{. 
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Figure 3: Quantile relative discrepancies: p = 3, q = 3, and power exponential distribution for Z{. 
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Table 4: Power of the modified tests when n = 20 and p = q = 2; normal, Student-t (i/ = 5) and power exponential (A = 0.6) 
distributions. 



normal Student-t power exponential 





7 = 


5% 


7 = 


10% 


7 = 


5% 


7 = 


10% 


7 = 


5% 


7 = 


10% 


V 


LK 


LRT 


LR* a 


LRT 


LR* a 


LRT 


LR* a 


LRT 


LR* a 


LRT 


LR* a 


LRT 


0.1 


5.4 


5.2 


10.6 


10.4 


4.8 


4.6 


10.1 


9.8 


5.0 


4.9 


9.8 


9.6 


0.2 


6.2 


6.1 


12.0 


11.8 


5.6 


5.5 


11.4 


11.2 


5.6 


5.6 


11.3 


11.1 


0.3 


9.5 


9.3 


16.4 


16.2 


8.1 


8.0 


14.8 


14.5 


8.3 


8.2 


15.1 


14.8 


0.4 


14.0 


13.8 


23.3 


23.0 


12.3 


12.1 


20.8 


20.6 


12.8 


12.6 


21.3 


21.0 


0.5 


22.3 


22.0 


33.4 


33.1 


18.4 


18.1 


28.7 


28.5 


19.3 


19.1 


29.9 


29.6 


0.6 


30.9 


30.6 


43.9 


43.5 


26.9 


26.6 


38.6 


38.3 


27.4 


27.1 


39.2 


38.9 


0.7 


43.0 


42.6 


56.5 


56.2 


35.3 


34.9 


48.6 


48.2 


37.0 


36.5 


50.0 


49.7 


0.8 


53.9 


53.6 


66.8 


66.4 


44.7 


44.4 


58.1 


57.8 


46.9 


46.5 


60.5 


60.2 


0.9 


64.9 


64.5 


76.0 


75.8 


57.3 


57.0 


69.8 


69.6 


56.5 


56.1 


69.3 


69.1 


1.0 


72.8 


72.5 


82.4 


82.2 


64.9 


64.5 


76.3 


76.0 


67.5 


67.2 


78.4 


78.1 


1.1 


81.5 


81.1 


88.4 


88.2 


71.1 


70.7 


81.0 


80.8 


74.5 


74.2 


84.0 


83.7 


1.2 


86.8 


86.7 


92.4 


92.2 


79.0 


78.8 


87.2 


87.1 


81.1 


80.9 


89.4 


89.2 


1.3 


90.9 


90.8 


95.4 


95.3 


84.4 


84.2 


91.0 


90.8 


86.0 


85.7 


91.9 


91.7 


1.4 


93.5 


93.4 


96.8 


96.8 


90.1 


89.9 


94.9 


94.9 


90.5 


90.3 


95.4 


95.3 


1.5 


97.4 


97.3 


99.6 


99.5 


93.2 


93.0 


96.9 


96.9 


92.3 


92.2 


96.4 


96.3 



5 Application 

The WHO MONICA Project was established in the early 1980s, and has been monitoring trends 
in cardiovascular diseases since then. One of the goals of this project is to relate changes in known 
risk factors (x) with the trends in cardiovascular mortality and coronary heart disease (y). In this 
paper, we consider a data set from the WHO MONICA Project that was considered in Kulathinal 
et al. [15], Patriota et al. [20], and Cao et al. [5]. This data set was first analyzed under normal 
distributions for the marginals of the random errors ( [15] , |20] ) , and thereafter under scale- mixture 
of normal distributions for the marginals of the random errors (Cao et al. [5]). The approach used 
in the present paper is different from the latter, since here we consider an elliptical joint distribution 
for the vector of the random errors. Although the class of scale-mixture of normal distributions is 
a special case of the the elliptical distributions, our proposal does not extend to the class proposed 
by Cao et al. [5]. This is because they assumed that the errors' distributions are independent, while 
we assume that they are uncorrelated but not independent. For our proposal, the errors will only be 
independent under normality. 

The data set considered here corresponds to the data collected for men (n = 36). As describe in 
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Table 5: Power of the modified tests when n = 40 and p = q = 2; normal, Student-t [v = 5) and power exponential (A = 0.6) 
distributions. 



normal Student-t power exponential 



V 


7 = 


5% 


7 = 


10% 


7 = 


5% 


7 = 


10% 


7 = 


5% 


7 = 


10% 


LR*. 


LRl* 


LR* a 


LR* a * 


LR* a 


LR* a * 


LRl 


LR* a * 


LRl 


LR* a * 


LR* a 


LRl* 


0.1 


7.4 


7.4 


13.6 


13.5 


6.5 


6.4 


12.4 


12.4 


6.7 


6.6 


12.6 


12.6 


0.2 


14.1 


14.0 


23.3 


23.3 


12.2 


12.1 


20.7 


20.6 


12.6 


12.5 


20.8 


20.7 


0.3 


26.5 


26.4 


38.2 


38.2 


22.5 


22.4 


33.1 


33.0 


23.0 


22.9 


33.7 


33.7 


0.4 


43.0 


42.9 


55.9 


55.9 


35.8 


35.8 


48.5 


48.4 


36.9 


36.8 


49.7 


49.7 


0.5 


60.7 


60.6 


72.4 


72.3 


51.4 


51.4 


64.4 


64.3 


53.6 


53.5 


65.9 


65.8 


0.6 


76.1 


76.0 


85.2 


85.1 


66.8 


66.7 


77.3 


77.3 


68.8 


68.7 


78.4 


78.3 


0.7 


87.5 


87.4 


92.9 


92.9 


79.1 


79.1 


86.9 


86.9 


80.3 


80.2 


87.8 


87.8 


0.8 


94.0 


93.9 


97.1 


97.1 


88.1 


88.1 


93.2 


93.2 


88.9 


88.9 


93.7 


93.7 


0.9 


97.6 


97.6 


98.8 


98.8 


93.7 


93.7 


96.8 


96.8 


94.0 


94.0 


97.1 


97.0 


1.0 


99.0 


99.0 


99.5 


99.5 


96.9 


96.9 


98.5 


98.5 


97.2 


97.1 


98.7 


98.7 


1.1 


99.6 


99.6 


99.9 


99.9 


98.5 


98.5 


99.3 


99.3 


98.8 


98.8 


99.5 


99.5 


1.2 


99.9 


99.9 


100.0 


100.0 


99.3 


99.3 


99.7 


99.7 


99.5 


99.5 


99.8 


99.8 


1.3 


100.0 


100.0 


100.0 


100.0 


99.7 


99.7 


99.9 


99.9 


99.8 


99.8 


99.9 


99.9 


1.4 


100.0 


100.0 


100.0 


100.0 


99.9 


99.9 


99.9 


99.9 


100.0 


100.0 


100.0 


100.0 


1.5 


100.0 


100.0 


100.0 


100.0 


99.9 


99.9 


100.0 


100.0 


100.0 


100.0 


100.0 


100.0 



Kulathinal et al. |15j . the data are trends of the annual change in event rate (y) and trends of the 
risk scores (x). The risk score is defined as a linear combination of smoking status, systolic blood 
pressure, body mass index, and total cholesterol level. A follow-up study using proportional hazards 
models was employed to derive its coefficients, and provides the observed risk score and its estimated 
variance. Therefore, the observed response variable, Y, is the average annual change in event rate 
(%) and the observed covariate, X, is the observed risk score (%). For these data, we use model @ 
with m = p = 1 and zero covariance between the errors e and u. 

Table [6] gives the two corrected statistics and the original likelihood ratio statistic for testing 
H.Q : (3 = 1 (p- values are given in parentheses). We considered the full sample (n = 36) and a 
randomly chosen sub-sample of n = 20 observations. We notice that for the full sample all the tests 
have similar p-values, and reject the null hypothesis at the 5% nominal level. For the sub-sample, 
however, the null hypothesis is rejected at the 10% nominal level when one uses either of the corrected 
statistics, but %$ is not rejected when the usual likelihood ratio test is employed. 
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Table 6: Corrected and uncorrected likelihood ratio statistics for testing T~Lq : (3 = 1; p-values are 
given in parentheses. 



Full sample Sub-sample 



LR 4.5037 (0.0338) 2.4160 (0.1201) 
LR* a 4.5701 (0.0325) 3.6583 (0.0558) 
LR** 4.5698 (0.0325) 3.5299 (0.0603) 



6 Concluding remarks 

We studied the issue of testing hypotheses in heteroskedastic multivariate regression models with 
measurement error under the elliptical class of distributions. For this class of distributions, the exact 
distribution of the likelihood ratio statistic under the null hypothesis is not easily attained and the 
usual asymptotic approximation to the chi-square distribution is required. However, for small sample 
sizes, this approximation may be very inaccurate and modified versions proposed by Skovgaard [23] 
can be used to obtain more accurate approximations under the null hypothesis. The main difficulty 
of this approach is in finding a suitable ancillary statistic. In this work, we find this ancillary statistic 
and all the quantities required to compute Skovgaard's adjustment. Simulation results indicate that 
the adjusted tests are more reliable than the original likelihood ratio test for small sample sizes. A 
real data application illustrates the usefulness of the results presented in this paper. 
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Appendix A. Observed information matrix 

The first-order derivative of Q with respect to 8j, with j = 1,2, ... ,s, is 



dl{8) 




i=l L 



(7) 



where W Po (u) — d\ogpo(u)/du . Define 




8 de/ 



1 rv 



a 



90/ 



da 



f; 



(i) _ d J^L 
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" 50, ' and _ 90, ■ 



Then, the vector fig. — dfi/dOj and matrix fig. = dVli/d9j are given, respectively, by 



Me, 



or 



0) 



and 



fi« 



T(3) 



0, 



0) 



for j, k = 1,2, ... ,s. Note that the derivative presented in ^ forms the score vector. The element (j, k) of 
the observed information matrix J = J(9) is Jg ,g k = —d 2 £(0)/d6jd9k, i.e., 



1 1 



2nJ k n~ 1 d l ~\ [djn^ng^di + 2nJ j ni 1 d^ +w P0 {dJ^7 1 d t )[dj(2n- 1 ng k ni 1 ng ] n- 1 



nr^Qg.g^-^di + 2/j,J k $l~ 1 Q,g j $X[ 1 di + 2 / iIfi7 1 fig fc fi-^ - fyle^T 1 * + 



for j,k — 1, 2, . . . , s, where Wp (u) = dW Po (u)/du. The vector fJ,g-g k — dfigJdOk is 

for j,k = 1,2, s. The matrix ^lg j g k — dflg j / dQk is symmetric and has dimension sxs, with j,k — 1,2, . . . , s, 
and is given by 



ttg 3 g k = c 



F (k) p TU) 



o 



p 0) pT(fe) 



for j,k — 1,2, ... ,s. 

In matrix notation, the observed information matrix is given by 



J=-T- R T M - V T ( 
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where the expressions for the elements of T, R, M, V, and Q are presented in Section [3] 

Appendix B. Derivatives with respect to the data 



We present the derivatives with respect to the sample space £', U' , and J. Consider the ancillary statistic 
a = (aj , aj , • • • , a^) T defined in (|6j), where dj = P^ 1 [zi — p,). Inserting Zi = Pidi + p, in the log-likelihood 
function Q, we obtain 

1(0; 9 , a) = ^ I -- log |fl<| + logp [(PiOj + £ - n) T Sir 1 (P iai +fi, — /j, 
Therefore, the j-th element of £' = <9£(0; 0, a) /90 is 

z=l 

Furthermore, the (j, fc)-th element of [/' = d 2 £(9; 9, a)/09d9 T is 

n ( T t 

I7j fc = 2^ I W P0 (dj^di) - (?e k ai + fa) n^fiT 1 * - fa.SlJ 1 (Pe k a % + fa 
i=i I L 

+ W' po (dJSl^di) (Pe^ + faY (-dj^ile^di - 2fan^d z ) I, 



Me* 



where P$. = dPi/9j, for j, k = 1, 2, . . . , s. We also have that the (j, k)-th element of J is 

~ n ( r t 

J jk = 2J2\W P0 (d[fl7^ - (Pe^ + fa) n^n^oj - j^.Qr 1 (p^a 
i=i I L 

+ (^fir 1 *) (p flfc a< + Jlg.yhr 1 ?^ (-aJP^h^ho^P^ - 2j5^Qr 1 P i a i ) j. 

In matrix notation, the derivatives with respect to the sample space used in the Skovgaard's adjustment are 

£' = 2R T w, U' = 2 (R T B + V T C) , and J = 2 (r t F + V T G^; 
with the elements of B, C, F, G, and w being defined in Section [3] 
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